C     POMSUB     Auxiliary subroutines (and functions).
C     Copyright (C) 1997  Filip Hroch, Masaryk University, Brno, CZ
C                    
C     This program is free software; you can redistribute it and/or modify
C     it under the terms of the GNU General Public License as published by
C     the Free Software Foundation; either version 2 of the License, or
C     (at your option) any later version.
C     
C     This program is distributed in the hope that it will be useful,
C     but WITHOUT ANY WARRANTY; without even the implied warranty of
C     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
C     GNU General Public License for more details.
C
C     You should have received a copy of the GNU General Public License  
C     along with this program; if not, write to the Free Software
C     Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
C
C    E-mail: hroch@physics.muni.cz
C    Web:    http://www.sci.muni.cz/~hroch/index.html
C    Postal address: F.Hroch
C                    Department of Theoretical Physics and Astrophysics
C                    Masaryk University
C                    Kotlarska 2
C                    611 37 Brno
C                    Czech Republic
C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

      SUBROUTINE SERAD(A)
C Sort array A from the greatest to the smallest.
      DOUBLE PRECISION A(3),T
      INTEGER I,J

      DO I=2,3
        DO J=3,I,-1
          IF( A(J-1).LT.A(J) )THEN
            T = A(J-1)
            A(J-1) = A(J)
            A(J) = T
          ENDIF
        ENDDO
      ENDDO
      END

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

      SUBROUTINE DATA(WATCH,SOUBOR,N,X,Y,IN,MMAX)
C Read data from file      
      CHARACTER*(*) SOUBOR
      INTEGER N,IN(*),NL,IP1
      DOUBLE PRECISION X(*),Y(*),P1
      LOGICAL WATCH
      
      OPEN(1,FILE=SOUBOR,STATUS='OLD',BLANK='ZERO',IOSTAT=N)
      IF( N.NE.0 )THEN
       PRINT *,' File not found :',SOUBOR(1:40) 
       N = 0
       RETURN
      ENDIF 
      CALL RDHEAD(1, NL, IP1, IP1, P1, P1, P1, P1, P1, P1, P1)

      N = 0
10    N = N+1
        IF( N.GT.MMAX ) GOTO 20
15      CALL RDSTAR(1,NL,IN(N),X(N),Y(N),P1,P1)
        IF (IN(N) .EQ. 0) GO TO 15
        IF (IN(N) .LT. 0) GO TO 20
      GOTO 10
20    CLOSE(1)
      N = N-1
      IF( WATCH ) PRINT *,N,' points was read from ',SOUBOR(1:40)
      END

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

      SUBROUTINE CENTER(N,X,Y,T)
C Calcualtion of mean arrays X,Y      
      INTEGER I,N
      DOUBLE PRECISION X(*),Y(*),T(2)
      
      T(1) = 0.0
      T(2) = 0.0
      DO I=1,N
        T(1) = T(1)+X(I)
        T(2) = T(2)+Y(I)
      ENDDO
      T(1) = T(1)/N
      T(2) = T(2)/N
      END

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC      

      FUNCTION ZOOM(N,X,T)
C                
      INTEGER I,N
      DOUBLE PRECISION ZOOM,X(*),T(2),ZV
      
      ZV = 0.0
      DO I=1,N
        ZV = (X(I)-T(1))**2
      ENDDO         
      ZOOM = ZV/(N-1)
      END
      
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC      

      DOUBLE PRECISION FUNCTION BMIN(N,X)
      
      INTEGER I,N
      DOUBLE PRECISION X(*),Y
      
      Y = X(1)
      DO I=2,N
        Y = MIN(Y,X(I))
      ENDDO
      BMIN = Y
      END
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC        

      DOUBLE PRECISION FUNCTION BMAX(N,X)
      
      INTEGER I,N
      DOUBLE PRECISION X(*),Y
      
      Y = X(1)
      DO I=2,N
        Y = MAX(Y,X(I))
      ENDDO
      BMAX = Y
      END

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

      SUBROUTINE USPOR(N,POLE,ID)
C Sort of array POLE by algorithm Quicksort
C By  N.Wirth : Algorithms + Data Structure = Programs 
C  N is dimension of array POLE and ID

      INTEGER N,ID(*)
      DOUBLE PRECISION POLE(*)

C This paramets give a size of stack. For M = 30 is a capacity 
C of array around half milions (5E5) elements in sorting array.      
      INTEGER M
      PARAMETER( M = 30 )
      
      INTEGER I,J,L,R,S,STACK(M,2),TMP
      DOUBLE PRECISION W,X
      
      S = 1
      STACK(1,1) = 1
      STACK(1,2) = N
10    L = STACK(S,1)
        R = STACK(S,2)
        S = S-1
20      I = L
          J = R
          X = POLE((L+R)/2)
30        IF( POLE(I).LT.X )THEN
              I = I+1
              GOTO 30
            ENDIF
40          IF( X.LT.POLE(J) )THEN
              J = J-1
              GOTO 40
            ENDIF      
            IF( I.LE.J )THEN
              W = POLE(I)
              POLE(I) = POLE(J)
              POLE(J) = W
              TMP = ID(I)
              ID(I) = ID(J)
              ID(J) = TMP
              I = I+1
              J = J-1
            ENDIF
          IF( I.LE.J ) GOTO 30
          IF( I.LT.R )THEN
            S = S+1
            STACK(S,1) = I
            STACK(S,2) = R
          ENDIF
          R = J
        IF( L.LT.R ) GOTO 20 
      IF( S.NE.0 ) GOTO 10
      END         

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

      SUBROUTINE USPOR2(N,POLE,ID)
C For DOUBLE PRECISION 'POLE'
C Sort of array POLE by algorithm Quicksort
C By  N.Wirth : Algorithms + Data Structure = Programs 
C  N is dimension of array POLE and ID

      INTEGER N,ID(*)
      DOUBLE PRECISION POLE(*)

C This paramets give a size of stack. For M = 30 is a capacity 
C of array around half milions (5E5) elements in sorting array.      
      INTEGER M
      PARAMETER( M = 30 )
      
      INTEGER I,J,L,R,S,STACK(M,2),TMP
      DOUBLE PRECISION W,X
      
      S = 1
      STACK(1,1) = 1
      STACK(1,2) = N
10    L = STACK(S,1)
        R = STACK(S,2)
        S = S-1
20      I = L
          J = R
          X = POLE((L+R)/2)
30        IF( POLE(I).LT.X )THEN
              I = I+1
              GOTO 30
            ENDIF
40          IF( X.LT.POLE(J) )THEN
              J = J-1
              GOTO 40
            ENDIF      
            IF( I.LE.J )THEN
              W = POLE(I)
              POLE(I) = POLE(J)
              POLE(J) = W
              TMP = ID(I)
              ID(I) = ID(J)
              ID(J) = TMP
              I = I+1
              J = J-1
            ENDIF
          IF( I.LE.J ) GOTO 30
          IF( I.LT.R )THEN
            S = S+1
            STACK(S,1) = I
            STACK(S,2) = R
          ENDIF
          R = J
        IF( L.LT.R ) GOTO 20 
      IF( S.NE.0 ) GOTO 10
      END         

    


CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

      SUBROUTINE ADD(NSTAR,ID1,ID2,NSTACK,STACK1,STACK2,SSTACK,S)
C Write next record or increased found serie
      INCLUDE 'munimat.inc'
                             
      INTEGER NSTAR,ID1(*),ID2(*)
      INTEGER NSTACK,STACK1(MAXSTACK,0:IDMAX),STACK2(MAXSTACK,0:IDMAX)
      INTEGER I,J,LL,JE1
      DOUBLE PRECISION SSTACK(MAXSTACK),S
      LOGICAL JE
                         
C Finding the serie in known series
      JE = .FALSE.
      DO I=1,NSTACK
        JE1 = 0 
        DO J=1,NSTAR    
          DO LL=1,NSTAR
            IF( ID1(LL).EQ.STACK1(I,J).AND.
     *          ID2(LL).EQ.STACK2(I,J) ) JE1 = JE1+1
          ENDDO
        ENDDO
        JE = JE1.GE.NSTAR-1
        IF( JE ) GOTO 70
      ENDDO 
70    IF( JE )THEN
C If found, increase its number 
        STACK1(I,0) = STACK1(I,0)+1
        STACK2(I,0) = STACK2(I,0)+1
        SSTACK(I) = MAX(SSTACK(I),S)
      ELSE  
C New record      
        NSTACK = NSTACK+1     
        IF( NSTACK.GT.MAXSTACK ) STOP 'Stack overflow.'
        STACK1(I,0) = 1
        STACK2(I,0) = 1
        DO J=1,NSTAR
          STACK1(I,J) = ID1(J)
          STACK2(I,J) = ID2(J)
        ENDDO
        SSTACK(I) = S 
      ENDIF  
        
      END

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

      INTEGER FUNCTION LN(A)
C Get position of last non-space character in A.
      CHARACTER*(*) A
      INTEGER I
      
      LN = 1
      DO I=1,LEN(A)
        IF( A(I:I).NE.' ' ) LN = I
      ENDDO

      END

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

      SUBROUTINE CTOF(A,B)
      CHARACTER*(*) A,B

      INTEGER I
      LOGICAL NULL

      NULL = .FALSE.
      DO I = 1,LEN(A)
        IF( .NOT.NULL .AND. A(I:I) .LT. ' ' ) NULL = .TRUE.
        IF( NULL )THEN
          B(I:I) = ' '
        ELSE
          B(I:I) = A(I:I)
        ENDIF
      ENDDO
      END

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

      SUBROUTINE GETPAR( MAXSTAR, NSTAR,CLIP, WATCH, BLANKA)
C Read parameters from file match.opt.                            

      INCLUDE 'munimat.inc'
      
      DOUBLE PRECISION X,CLIP
      INTEGER I,MAXSTAR,NSTAR
      LOGICAL WATCH,BLANKA
                    
      OPEN(1,FILE='match.opt',STATUS='OLD',IOSTAT=I)
      IF( I.NE.0 )THEN
        PRINT *,'File with parametres not found. Used default values.'  
        MAXSTAR = 13
        NSTAR = 5
        CLIP = 2.5
        WATCH = .FALSE.
        BLANKA = .TRUE.
        PRINT *,'Maximum read star = ',MAXSTAR
        PRINT *,'Number stars for identification = ',NSTAR
        PRINT *,'Sigma clipping factor = ',CLIP           
        PRINT *,'Watch progress = 0'
        PRINT *,'Blank line = 1'
      ELSE        
        READ(1,700,IOSTAT=I) X
        IF( I.NE.0.OR.X.GT.MMAX )THEN
          MAXSTAR = MMAX
          PRINT *,'Error during read parameters.',
     .            ' Used Maximum read star = ',MAXSTAR
        ELSE  
          MAXSTAR = X
        ENDIF
        READ(1,700,IOSTAT=I) X
        IF( I.NE.0.OR.X.GT.IDMAX )THEN
          NSTAR = 5
          PRINT *,'Error during read parameters.',
     .            ' Used Number stars for identification = ',NSTAR
        ELSE  
          NSTAR = X
        ENDIF
        READ(1,700,IOSTAT=I) X
        IF( I.NE.0 )THEN
          CLIP = 3.0
          PRINT *,'Error during read parameters.',
     .            ' Used Sigma clipping factor = ',CLIP
        ELSE  
          CLIP = X
        ENDIF  
        READ(1,700,IOSTAT=I) X
        IF( I.NE.0 )THEN
          WATCH = .FALSE.
          PRINT *,'Error during read parameters.',
     .            ' Used Watch progress = 0'
        ELSE  
          IF( X.GT.0.5 )THEN
            WATCH = .TRUE.
          ELSE
            WATCH = .FALSE.
          ENDIF    
        ENDIF
        READ(1,700,IOSTAT=I) X
        IF( I.NE.0 )THEN
          BLANKA = .TRUE.
          PRINT *,'Error during read parameters.',
     .            ' Used Blank lines = 1'
        ELSE  
          IF( X.GT.0.5 )THEN
            BLANKA = .TRUE.
          ELSE
            BLANKA = .FALSE.
          ENDIF    
        ENDIF
        CLOSE(1)
      ENDIF
      CLIP = CLIP**2
700   FORMAT(TR28,F10.0)      
      END    
        
           
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC    

